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We apply the Glasgow method for lattice QCD at finite chemical potential to a schematic random 
matrix model (RMM). In this method the zeros of the partition function are obtained by averaging 
the coefficients of its expansion in powers of the chemical potential. In this paper we investigate the 
phase structure by means of Glasgow averaging and demonstrate that the method converges to the 
correct analytically known result. We conclude that the statistics needed for complete convergence 
grows exponentially with the size of the system, in our case, the dimension of the Dirac matrix. The 
use of an unquenched ensemble at /i = does not give an improvement over a quenched ensemble. 

We elucidate the phenomenon of a faster convergence of certain zeros of the partition function. 
The imprecision affecting the coefficients of the polynomial in the chemical potential can be inter- 
^ peted as the appearance of a spurious phase. This phase dominates in the regions where the exact 

, partition function is exponentially small, introducing additional phase boundaries, and hiding part 

QQ ' of the true ones. The zeros along the surviving parts of the true boundaries remain unaffected. 
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CN ■ I. INTRODUCTION 

> 

In contrast with the numerous successes of lattice QCD, simulations at finite chemical potential |l|-|^Jl^ have 
- oscillated between mildly promising and outright frustrating. The source of the trouble lies in the following. In the 
QQ Euclidean formulation of QCD, the chemical potential spoils the anti-Hermiticity of the Dirac operator. As a result, 
. the fermion determinant is no longer a real number. In general, it has a complex phase. Hence, the action cannot 

serve as a statistical weight in a Monte Carlo sampling of field configurations. 
, Quenched simulations, where the fermion determinant is not included in the statistical weight, may provide a fairly 
reliable approximation to selected observables of the true unquenched theory. However, in the presence of a chemical 
potential, quenched simulations have produced consistently unphysical results |^. The reason is that the quenched 
theory is the Nf limit of an unphysical theory where the fermion determinant is replaced by its absolute value 
(jj ■ This is a theory with a second, "conjugate" set of anti-quark species together with the normal quarks. Because 

, of Goldstone bosons consisting of a quark and a conjugate anti-quark, the critical chemical potential in quenched 
' QCD is half the pion mass. This phenomenon was demonstrated in lattice simulations by Gocksch |^ using a U{1) 
. ^ toy model and was understood analytically in |^] by using a random matrix model inspired by QCD. 

One important conclusion of further studies of the same RMM is that the phase of the fermion determinant leads 
\ to very large cancellations in the ensemble averaging Q. A measure of this phenomenon is given by the fact that 
- - • the partition function is proportional to exp{—^^N), where N is the size of the random matrix, corresponding to 
the number of sites in a lattice simulation. Cancellations of this magnitude would require prohibitive statistics in 
order for a brute force simulatioij^ to be successful. We wish to note that in some models the construction of clever 
algorithms makes it possible to deal with these cancellations in Monte Carlo simulations |Q. 

It would be hard to overrate the potential importance of a successful lattice approach to QCD at finite chemical 
potential. At this time, there is not even an estimate for the value of the critical chemical potential from lattice 
simulations. In the absence of the guidance provided by the lattice, it is difficult to assess the many semi-empirical 
descriptions of nuclear matter at high density f^-p^] . For instance, an extension of the RMM to include temperature 
via the first Matsubara frequency Jisf gives reasonable predictions about the phase diagram of QCD, even while 
ignoring most of its dynamics. In the present paper we wish to exploit the qualitative similarity between our simple 
RMM and Nc = 3 QCD at finite chemical potential in an attempt to understand certain lattice results on the problem 
of finite /i. 



^ where the determinant is included as an observable in a quenched ensemble 
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We are interested in the analytic dependence of the QCD partition function on the chemical potential ^. This can 
be obtained by computing the coefficients of the expansion of the partition function in powers of the chemical potential 
or the fugacity. The Glasgow method of lattice QCD jlGj-pOt is designed to do this. The unquenched partition function 
can be seen as the quenched average of the fermion determinant, i.e. averaged only with the gauge action. In general, 
one may also use an unquenched ensemble at some fixed chemical potential fj, = ^q, and include the inverse fermion 
determinant at that same value, 

Zi,) = i^etPi,))„^c(^^) . (1) 

\ ^ \t^^> } I gauge, /^—/io 

Here C is an irrelevant constant. One expects the efficiency of the averaging process to depend on the overlap between 
the quantity being averaged and the distribution used to generate the ensemble. When these two functions have their 
largest values in vastly different places, this is known as an overlap problem. 

In the lattice Glasgow method the fermion determinant is expanded in powers of the fugacity ^ — exp(/i). The 
expansion is finite and exact, since the fermion determinant is just an N x N matrix {N is the number of lattice points 
times the number of degrees of freedom per site). It is obtained by writing det^ (/x) = det {P + £,), where P is 
called the propagator matrix. The expansion coefficients, written in terms of the eigenvalues of P, are then obtained 
by ensemble averaging. The zeros of the partition function in the complex fi plane map out the phase structure. In 
particular, the ones close to the real axis define the critical value(s) of fi. Since the original paper by Yang and Lee 
pl| this approach has been widely used in statistical mechanics. For some recent applications we refer to p^ ]. 

With the Glasgow method one obtains information about the full /u-dependence of Z, using an ensemble generated 
at one fixed value of the chemical potential. Of course, the question of the overlap between whatever ensemble one is 
using and the fermion determinant for the given /i remains. The issue is even more ominous considering that we are 
unable to even define the notion of an ensemble at nonzero /i, even in this random matrix model, due to the complex 
action. 

In this paper, we analyze the Glasgow method using a random matrix model at nonzero chemical potential. As we 
have already discussed before, this model mimics the problems of the QCD partition function at nonzero chemical 
potential. Since the phase structure of this model is known analytically, it is an ideal testing ground for evaluating 
this algorithm and obtaining a better understanding of its problems. 

In section 2 we introduce the random matrix model and derive some of its analytical properties. The bulk of this 
paper is in section 3. It contains our numerical analysis of the Glasgow method, and an explanation is given of the 
poor convergence of certain zeros of the partition function. Concluding remarks are made in section 4. 

II. RANDOM MATRIX MODEL 

We consider a random matrix model (RMM) defined by the partition function MS 



D{m,fj.) 



m iC + ^J. 
iC^ +11 m 



(2) 



Here, C is an x iV complex matrix, and the integration is over the Haar measure, 

/^^ = n / dQjdc:^ . (3) 

This model was first formulated [ p3|j2^ ] for /i = in order to describe the correlations of the smallest eigenvalues of 
the Dirac operator. In this case it has been shown rigorously that the model describes the zero-momentum sector of 
the low-energy effective partition function of QCD |^,^ . 

In the present context, the matrix D{m, fi) mimics the QCD Dirac operator for quark mass m and chemical potential 
fi. The integration over matrix elements replaces the integration over gauge field configurations. The massless part 
of our random Dirac operator is anti-Hermitian for fi — 0, but for fi ^ it has no definite Hermiticity properties. 
In order to study the properties of the partition function in the chemical potential plane /i, we rewrite the fermion 
determinant as follows, 
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detpK/i)) = det( + ^ ^^["+^ ) = dct(F(m)+/il) . (4) 



The matrix F{m) is analogous to the propagator matrix from lattice QCD [g7 28 1^ in the sense that its eigenvalues 
are the values of fi for which the fermion determinant vanishes. In terms of the eigenvalues Afc of F{m) the RMM 
partition function is 

Zim,fi) ^ J VCe-^'^^'^' l[{Xk + fi). (5) 

k 

The quantity (n) — In Z(m, /i) /N is the analog of the baryon number density of QCD. It is equal to the ensemble 
average of {j^ ^^^^^ ). In QCD at zero temperature, one expects the baryon number density to be identically zero 
for small fi, and then to increase starting from a certain critical value of fi p^ ]. Similarly, our model shows a phase 
transition with an increase of the baryon number density. However, (n) is not zero below the critical value of /i. See 
| p5[ for an explanation of the relationship between (n) and the baryon number density of QCD. 

It is not clear how to define a statistical ensemble of gauge field configurations (or random matrices for that matter), 
corresponding to the true partition function with finite fj,. In the quenched approximation one discards the fermion 
determinant, so the partition function does not depend any more on fj, or m. However, a 'number density' can still 
be computed by taking the average of X]fc xj+fj. ■ Similarly, one can define a quenched 'chiral condensate' from 
the ensemble average of -^Tr (Z?(to, . The quenched approximation can be interpreted either as the limit of a 
process where one takes the number of flavors to zero, or the quark mass m to infinity. The quantity 

k k ^ ^ 

is called the partially quenched baryon number density. The "sea" /xq defines the ensemble, and the "valence" /j, 
probes the eigenvalue distribution. 



A. Quenched eigenvalue distribution 

For m = the propagator matrix is block-diagonal, and its eigenvalues are i times those of C . The exact distri- 
bution of the eigenvalues of a general complex matrix has been calculated a long time ago by Ginibre . In our 
normalization, it is given by 

•' k>l 

The corresponding one-point function is 

P(A) =Ce--I^I^E^^ (8) 

fe=0 

In the large A'' limit the eigenvalues are uniformly distributed in the complex unit circle. This follows from the 
properties of the truncated exponential, to be discussed in more detail later. For |A| < 1, the truncated exponential 
is a good approximation of the complete one, so p is a constant. For |A| > 1, the truncated exponential behaves like 
a power of |Ap which is quickly suppressed by exp(— iV|Ap), so the distribution vanishes with a sharp tail of width of 
order 1/^/N. 

For nonzero m we can calculate the eigenvalue distribution of the propagator matrix in the large N limit using the 
conjugate replica trick We consider the partition function, where we have replaced the fermion determinant with 
its absolute value squared. 



det 



m tC + ^ 
iC^ + /i m 



(9) 



Here we can use either m or /i as a dummy variable probing the eigenvalue distribution of the corresponding operator 
as a function of the other variable which is made real and therefore has a physical meaning. For the spectrum of the 
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propagator matrix we set m = m* and probe the eigenvalues with the complex dummy variable ^. This is the reverse 
of what was done in |^ where the spectrum of the Dirac operator was investigated for given /i. The conjugate replica 
trick also allows us to calculate the eigenvalue distribution in the present case. 

Because of the absolute value of the determinant under the integral, the partition function is expected to be a smooth 
function of hf and the limit n/ — *■ should be obtainable from the partition function for positive integral values of 
fif. For any fif, the eigenvalue density is positive definite so it cannot be an analytic function of a complex variable. 
The resolvent is defined as G{z) = jji^f. ^ "L^ )- The average eigenvalue density is then given by NRe {dz'G(z)). For 
the partition function the resolvent in the fi plane is G{z = —fi) — d^Z{m, fi, ^*). 

By standard manipulations |^ the partition function can be rewritten as an integral over an 2ny x 2hf complex 
matrix. For hf ~ 1, one finds 



Z{m, m*, fi, ^*) 



(fa d% (fc <Pd e-^('^'^*+''^'*+-*+'i'i*) det 



" m + a 







id 




m + a* 


ic 








id* 


m* + b* 




ic* 







m* +b 



(10) 



The resulting saddle point equations have two kinds of nontrivial solutions, depending on whether the off-diagonal 
quantities c, c*, d, d* are identically zero or not. If they vanish, the partition function factorizes into pieces that depend 
only on m, or on w*, Then, G{z) is always an analytic function of z. Therefore, the region in the complex plane 
of z (m or fi) where the eigenvalues are located must be dominated by the other kind of solution, which mixes the 
parameters and their conjugates. It turns out that this solution is such that c = c* = d ^ d*, and the quantity cc* is 
positive. For m = m* real, it is given by 



7 ? 



(11) 



the 



The boundary of the domain of eigenvalues is given by the curve, where the two types of solutions merge, i.e 
condition cc* = 0. One can see immediately that for m = the boundary cc* = reduces to the unit circle. 

In Fig. |l] we show the distribution of the eigenvalues in the complex plane of the propagator matrix of size A'^ = 96 
for masses of value m — Q and m = 0.0625. Each plot consists of eigenvalues from 20 independent configurations. 
The analytic curve given by cc* = is also drawn in each case. We clearly see that this does give the correct result 
for the boundary of eigenvalues. 




'S: 



-1 -0.5 0.5 1 -1 -0.5 0.5 1 



(a) m = (b) m = 0.0625 

FIG. 1. Distribution of the eigenvalues in the complex plane of the propagator matrix of size A'' = 96 along with the analytic 
curve for the boundary given by cc* = in (0). 
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B. Unquenched partition function 



The unquenched partition function defined in eq. (||) can be computed analytically It is given by 

Z{m, ^i)^ f dada^e-''''^^' det f + ^ (12) 
J y fi m + J ^ ' 

where the integration is over the Nf x Nf matrix a. In the N oo limit the integrals can be evaluated via saddle 
point approximation. The saddle points are given by a cubic equation, 

fT = (7 

a{m + aY — /i^cr = m + cr , (13) 

where the matrix variable is diagonal and cr is now a number. For fixed real m, there are four branch points in the 
complex [I plane, given by four of the six [i, values for which the discriminant of the above equation, 

D - (™'^' (^^' " " i) + + ^')') ' (^^^ 

vanishes. The branch points are connected by two branch cuts. The derivatives of the partition function are dis- 
continuous across these cuts. The points where the cuts cross the real axis are interpreted as critical values of the 
parameters. They indicate a first-order phase transition in the thermodynamic limit pj6[|r^. 

For finite A^, the partition function can be evaluated as a polynomial in either m or /i. The coefficients arc obtained 
by expanding the determinant and performing the integrals. One exact expansion of the partition function is given 

by [| 

^("^'A^) - ]^l.Z.^fc!)2 {N-3~k)\ ■ ^^^^ 

The zeros of this polynomial can be readily calculated. They are located along lines in the complex plane and in the 
limit — > oo they converge to the exact cuts found by the saddle point analysis below (^. 

C. Phase structure and zeros of a partition function 

When a partition function has a nontrivial phase structure in the thermodynamic limit, the complex plane of one 
of its thermodynamic parameters is split into regions separated by cuts. Inside each region, the partition function is 
analytic in the parameters, so that its derivatives are continuous. The (first-order) transitions occur across the cuts, 
where the (logarithmic) derivatives of the partition function will in general have discontinuities. 

Before taking the thermodynamic limit (i.e. at finite iV), the partition function is an analytic function of the 
parameters. If this partition function at finite is a polynomial in any of its thermodynamic parameters, such as in 
our case the mass or the chemical potential, it will have zeros in the complex plane of these parameters. Its logarithmic 
derivatives will have poles at the locations of the zeros. In the thermodynamic limit, the zeros coalesce into the same 
cuts that define the phase structure. 

Our partition function illustrates nicely the connection between the zeros and the cuts. For = it is a truncated 
exponential, 

^HVy .... 

For /i^ ^ 1 the largest term in this sum occurs well before the truncation so we have a good approximation of the 
exponential. On the other hand, for /i^ ^ 1 the series is dominated by the term with j = N . One can obtain a better 
estimate using the incomplete gamma function, which is closely related to the truncated exponential [ pl| : 

oo 

n ^ /• 

n!e-"^||- = V{\ + n,x) = e-*edt . (17) 
k=o ' 
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For real positive arguments the situation is simple. If x is less than the value t = n which maximizes the integrand, then 
the saddle point is integrated over, the integral is a good approximation to the gamma function and the exponential 
is obtained. For x > n the integral is dominated by its endpoint and the integrand is well approximated by the last 
term of the series. 

If X is complex, the saddle point dominates if the integration contour can be deformed across the saddle point. For 
our partition function 



'N 



du , 



(18) 



this happens if Re(/i^ + In/i^) < —1 and Re/z^ > —1. In this case ZAr(0,/i) ^ exp(— /i^TV — N). If the integral is 
dominated by the endpoint then ZAr(0, /i) ~ /x^^. Nontrivial zeros of the partition function are obtained for values of 
the chemical potential where the saddle point contribution and the end point contribution arc of comparable order of 
magnitude. The zeros are thus given by the equation [|| 



Re(l + Ai^ + log/i2) = 0. 



(19) 



At finite N the deviation of the zeros from this critical curve is of order 1/vN 
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FIG. 2. Zeros from averaging over 10^ and 10^ configurations and exact zeros for a quenched ensemble with — 32. 



III. GLASGOW AVERAGING IN THE RMM 



The unquenched partition function for given /x and m can be thought of as the quenched expectation value of the 
fermion determinant. If the zeros of the partition function are known, we can express the partition function as a 
product over its zeros 

(U_(>^k + ^^)) = ll(^k + ^^) ■ (20) 

k k 

In this identity, the zeros of the partition function appear as a kind of "averaged eigenvalues" of the propagator 
matrix. The Glasgow method from lattice QCD attempts to perform this averaging. For a given matrix from the 
quenched ensemble, one writes out the coefficients of the polynomial on the left-hand side. The average coefficients 
are the coefficients of the right-hand side. The main focus of this paper is to study the convergence properties of the 
zeros to the exact ones for finite ensembles. 

The eigenvalues are in general complex. Because of the structure of F(rn), the eigenvalues occur in pairs {A^, — A^}. 
Also the matrix — C occurs with the same probability as C in (^). Therefore the ensemble average will also contain 
the pair {— Afe,A^}. Upon ensemble averaging, one can then easily show that the odd coefficients vanish and the 
remaining coefficients are real. 
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These simplifications should not mislead us into believing that we have safely avoided the trouble of averaging 
over the complex phase of the determinant. As we will see shortly, the problem will show up in the form of a very 
high precision needed for the coefficients in order to calculate the roots reliably. Since the suppression achieved by 
averaging over the phase of the determinant is on the order of the magnitude of the unquenched partition function, 
which for m = is exp (— iV(l + /i^)), it becomes exponentially difficult to achieve such precision. 

The quenched ensemble is not the only way to sample the set of all matrices. One may multiply and divide by any 
convenient function in the above formula. One factor modifies the ensemble, i.e., the way the individual matrices are 
generated, and the other one is used to compensate each configuration for the modified weight. One obvious choice 
is to use the unquenched ensemble at ii — 0. 

In the next section we report the results of numerical experiments performed using Glasgow averaging for the 
partition function (0), both with the quenched ensemble and the unquenched ensemble at ^ — 0. 
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3. Zeros from averaging over 10^, 10® and 10* configurations and exact zeros for ensembles with A'^ = 16. 



A. Numerical simulations 



We performed two different kinds of simulations. For m 7^ 0, we generated matrix elements with the Gaussian 
distribution given by (^). We constructed the propagator matrix and obtained its eigenvalues. For each set, the 
coefficients of the corresponding polynomial are calculated and added to the average. For m = 0, a much more 
economical procedure is possible, namely, generating sets of eigenvalues directly, using the exact Ginibre distribution 
. In this case we have employed a Metropolis algorithm. A given set is varied using small steps and the modified 
set is accepted depending on the corresponding value of the weight function. We have found that both cases have 
similar convergence properties, and in this paper we will only report on the case to = 0. 

In Fig. we show the zeros in the complex /x plane for a quenched ensemble of matrices of size = 32 averaged over 
A'^; = 10 ,10'' configurations. The averaged roots do converge to the exact values as expected, but the convergence is 
extremely slow. The unconverged zeros are the ones situated closer to the real axis. They are situated in a cloud of a 
well defined shape and most of them are located on its edge. As the averaging proceeds the cloud shrinks and finally 
disappears. All the roots situated outside the cloud are obtained correctly for a given A^, A^^ combination. This is 
illustrated in Fig. |^ where we perform ensemble averaging as large as Ne = 10®. In order to obtain more converged 
roots, we had to reduce the size of the matrices to A = 16. In this figure we also show results for Ne = 10^ and 
10^. The roots next to the real axis are always the last to converge. This is unfortunate since the value of the critical 
chemical potential (in our case fic = 0.527) is determined precisely by discontinuity across the real axis. The same 
pattern is observed for various matrix sizes A^. The shape of the cloud of unconverged zeros is very similar. However, 
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the number of configurations Ne corresponding to a given degree of convergence increases sharply with the matrix 
size N. One is able to determine fic with reasonable accuracy only for small values of N. As we can see in Fig. H, 
even in that case a very large number of configurations is required. Only after averaging over 10^ configurations do 
we get a good estimate of pic, but the roots are still not completely converged. 




FIG. 4. The width A{Ne) of the unconverged cloud as FIG. 5. The logarithm of the number of configurations 

a function of the number of configurations Ne- Notice (Ne) needed to achieve the given error {A{Ne)) in the 

that A{Ne) shows a logarithmic dependence on A''^;. zeros near the real axis for different sizes of matrices (A*'). 



Convergence is not improved by including the determinant at ^ = in the statistical weight used to generate the 
eigenvalues. It was suggested previously Q that using an unquenched ensemble at zero chemical potential might 
improve the efficiency of the averaging. Our results indicate that, at best, using the unquenched ensemble has no 
effect, but as it appears from Fig. ^the results are actually worse than for the quenched ensemble. An easy explanation 
is that a small number of configurations with the determinant close to zero are assigned a very large weight factor 
(cf. (|^)) which has a destructive effect. Using a "negative" number of fiavors is another possibility but it is not likely 
to improve convergence. 

We were able to run simulations which show some noticeable degree of convergence with N up to 48. There are 
many ways one could measure the degree of convergence of a given set of roots. The measure we will employ is the 
distance A{Ne) between where the two curves forming the boundary of the unconverged roots cross the real axis. 
For the lower boundary we use the smallest imaginary part of the eigenvalue near the real axis. The upper estimate 
is made by running a circle through the points near the real axis and choosing the one with the largest radius. In 
both cases the determination of how close to the real axis the points should be is made by eye, choosing a cutoff that 
appears to give a reasonable fit to the boundary near the real axis. Ideally these are the points where the contour of 
the cloud of unconverged roots crosses the real axis, corresponding to a lower and an upper estimate for fic using the 
given ensemble. 

The dependence A{Ne) on the size of the ensemble, Ne, is illustrated in Fig. |^ where we show results for N — 12, 
TV = 24 and N = 48. The thickness of the cloud of zeros, A{Ne), shows a logarithmic dependence on Ne- We have 
included logarithmic fits for the range Ne = 10^ through 10® (10^ for N = 12). Of course, once A{Ne) becomes close 
to zero, it stays that way. The slopes of the fits for N = 24, 48 vary only slightly with N- 

Finally, the issue of most practical interest is how the number of configurations Ne required to achieve a fixed 
value of A{Ne) varies with the size of the matrix. This was estimated by fitting lines to the available data for A{Ne) 
versus logig(iV£;). As we saw in Fig. |4|this does give nice linear fits. We then chose values of A(7V£;) and calculated 
where the linear fits intersected these values for each value of N- Our results are plotted in Fig. H. We conclude 
that the number of configuration required to obtain a given precision increases exponentially with the matrix size 
N. This is the most important conclusion of this paper. In the remaining sections we will attempt to reinforce this 
conjecture by studying the sensitivity of the zeros of the exact partition function to small random perturbations of 
the corresponding polynomial coefficients. 
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B. Perturbed exact zeros 



It is a well known fact that extreme care has to be taken when calculating zeros of very high order polynomials. 
In particular, this is the case for a finite representation of a partition function, such as the one under investigation, 
where one is ultimately interested in large values of N. In our numerical experiments the polynomial coefficients are 
obtained as a result of a statistical averaging process. The averaged coefficients are approximations of the exact ones. 
As the simulation proceeds, the error decreases. 

We wish to investigate the sensitivity of the zeros of the partition function to the precision with which the coefficients 
are obtained. We consider the exact polynomial and add a fixed relative error to each coefficient: 



Ck = Cfc(I + Rke) 



(21) 



where Rk are random numbers between [—1, 1], and e is a small real positive number. We then calculate the zeros of 
the polynomial ^ c^-z'''. 
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FIG. 6. Zeros from averaging over Ne = 10^ configurations and zeros of the exact polynomial perturbed with an error 
factor of e = 10~^ both for A'' = 32. The pattern of the false zeros is practically the same for Glasgow averaging and the 
polynomial with artificial noise. 



The effect of obtaining approximate coefficients by Glasgow averaging is quite similar to that of simply adding 
"artificial noise" to the exact coefficients. In Fig. |6| we plot the roots obtained for = 32 via Glasgow from Ne = 10^ 
configurations and the roots of the exact polynomial perturbed by noise of magnitude e = 10"'^. The patterns of the 
two sets of roots are hardly distinguishable. 

In this subsection we study the effect of such small random perturbations. In particular, we are interested in the 
dependence of the precision A{Ne) (in obtaining fi^) on TV and e. This approach has the advantage that we can 
consider larger values of N than in the case of direct Glasgow averaging. We expect that the relation between the 
error parameter e and the number of Glasgow configurations necessary to achieve the same accuracy is given by the 
central limit theorem, e ^ 1/ \/Ne- 

One striking feature for larger N is the extreme sensitivity of the roots to small perturbations. For N — 96, a noise 
factor of 10~^* already leads to a 20% error in the critical value of the chemical potential. Therefore, computing the 
zeros for = 96 is already beyond the capability of standard double precision.^ 



^ In all our calculations we use multiprecision arithmetic, implemented either using the GNU multiprecision package or the 
one made publicly available by NASA [p3|. 
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In Fig. we plot zeros with different degrees of artificial noise (see label of figure) for N — 96. The zeros for 
e — 10~^^ coincide with the exact roots. The noise factor for the remaining four sets of zeros increases with the 
same factor of 10"^ from one set to the next. The intercepts of the edge of the eigenvalue cloud with the real axis are 
practically equally spaced for the different sets of zeros. The precision A(e) appears to be a linear function of log(e) 
until A(e) vanishes. 

In Fig. ||we plot A(e) versus — log]^Q(e) for N — 48, 96, 192, 384. For error values not too for away from convergence, 
we observe a linear relationship between A(e) and — log(e). In all cases convergence is approximately achieved when 
— \ogiQ{e) « N/4:. This shows that the accuracy required to achieve a given precision increases exponentially with 
N. These results are consistent with results obtained with Glasgow averaging where the relative variance of the 
coefficients is roughly constant and varies as ~ I/^/N e with only a weak dependence on N . 

We conclude that adding random noise to the exactly known coefficients has the same effect on the roots of the 
partition function as Glasgow averaging. This confirms once more that an exponentially large number of configurations 
is required to obtain a given precision. Extrapolating the A'^-dependence of the number of configurations needed given 
by the linear interpolation in Fig. |^ to lattice simulations leads to extremely large numbers for even a small lattice 
size. Consider for example a matrix size of = 128 (so D in (||) is 256 x 256), corresponding to a 4"^ lattice with one 
degree of freedom per site. To achieve a precision of A = 0.3 we would need approximately Ne — 10^^ configurations, 
and for A = 0.2 we would need Ne = 10^^ configurations. A much more reasonable number of configurations can only 
be achieved if we consider a small lattice of 2^ which corresponds to = 24 for Kogut-Susskind (staggered) fermions. 
Here we would need Ne = 10^ configurations for A = 0.3 or Ne = 10^ configurations for A = 0.2. This exponential 
dependence was noted previously 011 • In 0' 

a signal was achieved for a lattice of 2"* with 10^ configurations but 
for a 4"^ lattice only a very weak signal was found after 2 x 10^ configurations. 



C. Further analysis of the convergence of Glasgow averaging 



In the previous two subsections we hope to have convinced the reader that the phenomena accompanying Glasgow 
averaging are nothing but the effect of knowing the polynomial coefficients of the partition function only with limited 
precision. This follows from the fact that the pattern of the false (unconverged) zeros is reproduced by adding small 
random numbers to the exact polynomial coefficients. 

In our model the zeros close to the real axis are the last to converge. These are also the most interesting in practice, 
since they determine the critical value /j,c. It is conceivable that for certain situations in QCD (such as with finite 
temperature) the roots close to the real axis are among the faster converging ones, which would provide a glimmer of 
hope for the Glasgow method. 



10 



In this subsection we consider the effect of artificial noise in detail. We wish to understand this phenomenon as 
well as why the false zeros are always concentrated in a cloud of a well defined shape. We will find that there is a 
clear correlation between the magnitude of the partition function and the stability of the zeros. 




-1 -0.5 6 0.5 1 

FIG. 9. Level curves of the absolute values of the free energy per site for the RMM as a function of the complex parameter 

fi. 

In Fig. |9| we give a topographic map of the absolute value of the RMM partition function on a logarithmic scale as 
a function of fi. That is, we plot the level curves of the free energy per site given (up to a constant) by ln(|ZArp)/7V. 
This is the natural quantity to study, since the partition function scales exponentially with N |^,^. We observe a 
discontinuity in the derivative of the absolute value of the partition function in the complex fj, plane. This is the locus 
of the zeros of the partition function. The deepest points are those around the (real) critical value of the chemical 
potential. A quick comparison with the preceding scatter plots should convince the reader that the higher the value 
of In(l^l) the more robust are the zeros in that neighborhood. However, the level curves do not match the contour of 
the cloud of false zeros. The situation a is little more subtle as we will see below. 



1. The noisy partition function 

It is useful to separate the approximate ( "noisy" ) partition function into the exact one and the part due to the 
error in the coefficients, 

Ztot{y) = Zo(^) + Zcrr(M) ■ (22) 

The partition function Zcrr(M) is a polynomial whose coefficients are the differences between the exact coefficients and 
the approximate ones, 

N 

■^crr(Ai) = '^Sck^Ji^ ; Ck^Ck+6ck . (23) 
fc=0 

The differences are generated as 5ck = ^RkC-k in the 'artificial noise' case. Our initial assumption was that for true 
Glasgow averaging the relative error in the coefficients is approximately the same for all zeros. This assumption was 
reinforced by the similarity between the two types of results discussed in previous subsections. In the following, we 
will discuss mainly the 'artificial noise' partition function. 

In terms of Zorr(A')) an explanation of the qualitative picture we have observed is the following. In the regions of 
the complex /i plane where |Zorr(M)| ^ I^oIa*)! one may ignore the Zd-T^fJ-). The roots of the total partition function 
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located in these regions coincide to a good degree with the corresponding exact roots. On the contrary, in regions 
where Zcr^ifJ-) dominates, the roots are determined by the latter and their general pattern has no similarity to their 
exact counterparts. This is consistent with the fact that the robust zeros are the ones situated in the region with 
higher ln(|Zo|). The shrinking of the cloud of false zeros with decreasing e is a consequence of the corresponding 
decrease in the magnitude of Zcrr- 

Throughout most of the complex /i plane, Zq is clearly larger than Zcm-- However, near the roots of Zq the error 
part has its chance to dominate. This is because the value of the exact partition function is the result of a major 
cancellation in the region where Zq is well approximated by exp(— ^^A^). The sum of the (finite) series is much smaller 
than the general term, which is typically of order 1. By multiplying each term of this series with a random number 
we spoil the cancellations. The sum of the perturbed series may therefore be significantly larger in absolute value 
than the exact sum. 




2. Sensitivity of the zeros 

One arrives at similar conclusions by studying the infinitesimal variation of the roots. Let fik be the exact roots 
(so that Zo(/iA;) = 0) and let /ifc + SfXk be the roots of the total noisy partition function. From our decomposition, we 
have 

Zoi^^, + s^^,) + z^U^^, + s^,,) = - (24) 

Since the modified roots are not zeros of Zq or Z^m the above equality is not fulfilled trivially. It indicates that the 
false zeros should be located in the transition region where \Zq\ and [.Z^errl are comparable, i.e., on the border of the 
region where Z^rr dominates, rather than scattered inside it. In Fig. we show how the contour of the cloud of false 
zeros is obtained using the formula above. 

We can also make a quantitative estimate of how well converged a given root is. If we expand the first half of ( |2^ ) 
we get 

Zo{^J-k) + SUkZoipk) + 2'orr(Mfc) + '5A*fc^orr(Mfc) ~ 0. (25) 

The first term is simply zero. The last term we will neglect as being of higher order in dfik- The remaining two terms 
can be rearranged to give 
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S,^, « - (26) 

In other words, the variation of a given root fik is proportional to the ratio of Z^rr to Z'q. If this quantity is neghgible, 
the root is close to the exact one. If the ratio is close to 1 or larger, the shift in fik is large, and the root is not 
obtained correctly. Therefore we can say that in the region where l^orrl l-^ol' roots of Ztot are reliable. 




3. Error piece as additional phase 



We may check our arguments in the previous sections by comparing the magnitudes of the exact, the total and the 



noise partition functions (or free energies) in the complex n plane. Fig. 11 shows the absolute value of the free energy 
per site, hi(Zo)/N in one quadrant of the complex /i plane. The location of the exact zeros is on the cusp which runs 
from fi = —i to fj, = 0.527. 

In Fig. |l2| we have the same type of plot now for the total partition function with a noise factor e — 10^^. The cusp 
in the noisy result bifurcates. The exact and the noisy surfaces coincide exactly except for the region between the two 
new branches of the cusp, where the noisy partition function is larger. The new cusps coincide with the locus of the 
false zeros. A few false zeros are also scattered inside the region. It is clear from Fig. ^ where we plot the surfaces 
corresponding to Zq and Zorr simultaneously that the new cusps are located at the intersection of the free energy 
surfaces. The exact piece dominates everywhere except for inside the bifurcation, where the error piece dominates. 
The /i dependence is so steep for both of them that the smaller piece becomes negligible very fast as one moves away 
from the intersection line. 

When we discussed earlier the properties of polynomials as partition functions, we associated the different analytic 
functions which are approximated by the polynomial in different regions of the complex parameter space with the 
phases of the partition function. We mentioned that the zeros are typically located along the phase transition lines, 
since within any given phase the partition function is a smooth analytic function which does not vanish. The partition 
function Z^rr has a quasi-analytic behavior similar to that of the partition function itself. The bifurcation of the line 
of zeros can be interpreted as the presence of an additional spurious 'phase' in the partition function, namely, the 
region where the error piece dominates. Then the fact that the false zeros are located on the \Zcri\ = \Zq\ line is 
natural. 
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4. Analytic approximation of the error partition function 

The question is whether the error partition function can be at all approximated by an analytic function. In the 
present case the answer is very simple. Let us consider the error partition function, 



N 



(a.) = ^E^fc ' ■ ^^^^ 

fe=0 

It is the same truncated inverse exponential we have discussed, only now each term is multiplied by a random number 
Rk of order 1. The terms in the original series conspire to achieve a major cancellation, of order e~^^^ . The random 
numbers spoil this, and each term in the series is left to fend for itself, and the sum is dominated by the largest terms. 
From the ratio of two consecutive terms in the sum, 

we conclude that the largest term (in absolute value) is the one with k = kj^ax — l/^^^l- If Im^I > 1 then the largest 
term is the one with k — N . Of course the surrounding terms must have a significant contribution, but the end result 
should still be proportional to this largest term. For < 1 our estimate is therefore 

« ^ -—= (29) 



In Fig. 1^ we plot the absolute value of the two surfaces corresponding to the continuum hmit, e^''^^ and e^/i^^, 
and the one corresponding to our estimate of the error partition function given above. The intersections of the three 
surfaces follow the pattern of the corresponding zeros. We also checked that (|29| ) approximates the error partition 
function well. 



14 



4 




FIG. 14. The exact and the error partition functions computed analytically. 



As an added bonus, we can now explain the scaling with N of the required precision. We found that the error 
partition function, i.e., the exact partition function whose coefficients have been multiplied by random numbers of 
order 1, is well approximated by an analytic expression. This expression shares an important property of the true 
partition function, namely, it scales exponentially with N. 

i-ln(Z— *(^)) « -lln(Z— ; i-ln(Z-(A.)) « i-ln(Z-(A.)) (30) 

The magnitude of the error partition function is also controlled by the factor e which mimics the precision to which 
the coefficients are approximated in the averaging process. The locus of the false zeros is controlled by the relative 
magnitude of the true partition function and the error partition function. To have the same pattern of zeros, we must 
also scale e exponentially with N, e{N) = . 



IV. CONCLUSIONS 



We have investigated Glasgow averaging using a random matrix model at nonzero chemical potential. We have 
found that in a quenched ensemble the method converges, but that it requires an exponentially large number of 
configurations. 

The roots of the averaged polynomial are initially distributed similarly to the eigenvalues of individual configu- 
rations. As the averaging proceeds, the roots approach their exact values. After averaging over a finite number of 
configurations, the roots clearly separate into two groups. Some roots are close to the corresponding exact ones. Typ- 
ically, these are the zeros far from the real axis. The remaining roots are situated in a cloud around the intersection of 
the real axis and the locus of the exact zeros. The zeros outside the cloud are practically exact, while those inside and 
on the boundary of the cloud are very badly determined. They cannot be traced to individual exact zeros. The shape 
of the cloud is similar for different matrix sizes. It shrinks as more configurations are taken into account. However it 
only shrinks as a logarithmic function of the total number of configurations. 

By interpolating the number of configurations needed to reach a given precision for several matrix sizes, we were 
able to estimate the dependence of the number of configurations needed for a given matrix size. Our conclusion is 
that the number of configurations needed to reach a given level of precision grows exponentially with the size of the 
matrix. 

The results of the corresponding unquenched simulations, using an ensemble generated with Nf = 1 and fi = are 
similar to the quenched ones. The unquenched ensemble may have been helpful in improving the overlap between the 



15 



simulation and the 'true' ensemble corresponding to fixed nonzero /x. However, this would still necessitate simulating 
an ensemble at nonzero /x, which is precisely what we were trying to avoid. The exponential statistics observed by us 
are more likely to be the signature of the sign problem itself, i.e., the magnitude of the cancellation brought about by 
the varying phase of the determinant. Hence the Glasgow method is unable to surmount the sign problem. However, 
in a problem where the latter is absent, such as SU{2) simulations |34|, the Glasgow method - even quenched - should 
be helpful. It would be interesting to see how the overlap issue manifests itself in this case. 

We obtained results very similar to those of Glasgow averaging by perturbing the exact coefficients in the RMM 
partition function polynomial. We studied empirically the dependence of the reliability of the zeros on the precision 
with which the coefficients are known and on the size N of the model matrix for larger matrices. Just like in the case 
of averaging, this dependence is exponential. That is, given a desired error bound for the zeros, the necessary precision 
on the coefficients grows exponentially with TV. The extrapolation of this dependence to Glasgow averaging translates 
into exponentially large statistics, since the precision on the coefficients should be proportional to the square root of 
the number of configurations. 

There is a correlation between the phenomenon of slower or faster converging zeros and the magnitude of the 
continuum partition function. The zeros that converge slowly are in the region where the partition function is 
suppressed. Large cancellations require better precision, hence more statistics. More formally, the effect of perturbing 
the coefficients can be understood as the addition of an extra (error) polynomial to the true partition function. This 
extra piece is found to scale exponentially with N, just like the true partition function. The effect on the phase 
structure can be seen as the introduction of a spurious phase, which replaces the true ones in the regions of parameter 
space where the error partition function dominates. The zeros in these regions follow the modified phase boundaries. 
The scaling property of the error partition function explains the need for exponential statistics in Glasgow averaging. 



The negative result regarding unquenched simulations at /z = is perhaps disappointing. It indicates that the 
quenched ensemble and the unquenched ensemble at /it = are equally relevant to the behavior of the model close to 
the critical /it. The upside of the failure of the unquenched ensemble at /i = is that in future applications of the 
Glasgow method it might be worth trying to use a quenched ensemble. We remind the reader that using the quenched 
ensemble is not equivalent to the quenched approximation. 

Another glimmer of hope is in the phenomenon of fast converging zeros. In the present random matrix model 
the zeros that determine the critical /i are the last to converge. From our analysis there is no indication that the 
sensitive zeros are generally those close to the real (physical) axis. There is a possibility that in QCD or in other 
interesting non-Hermitian models the critical parameter values are determined by the robust zeros. It is hard to make 
any statement in this respect from the currently available QCD data p6[ |. Perhaps a schematic model with features 
closer to QCD in terms of eigenvalue distribution in the complex /i plane would help clarify this issue. In general, it 
would be interesting to see how much of the analysis in this work regarding zeros of approximately known polynomial 
partition functions applies to other models. 
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